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We perform 3+1 general relativistic simulations of rotating core collapse in the context of the 
collapsar model for long gamma-ray bursts. We employ a realistic progenitor, rotation based on 
results of stellar evolution calculations, and a simplified equation of state. Our simulations track self- 
consistent ly collapse, bounce, the postbounce phase, black hole formation, and the subsequent early 
hyperaccretion phase. We extract gravitational waves from the spacetime curvature and identify a 
unique gravitational wave signature associated with the early phase of collapsar formation. 
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There is strong observational evidence linking long 
gamma-ray bursts (LGRBs) with the death of massive 
stars in core collapse (e.g., pQ). It appears likely that 
LGRBs are made in metal-poor progenitors with de- 
generate iron cores. These may be ordinary massive 
stars turned into Wolf-Rayet objects by mass loss or bi- 
nary interactions [2, 3 or, perhaps, peculiar, fully- mixed 
stars [U [5] . Both could result in a type-Ibc core-collapse 
supernova (CCSN) harboring a LGRB central engine. 
The nature of the latter and the details of the CCSN- 
LGRB relationship are uncertain. Viable engine set- 
tings all require rapid progenitor rotation and include 
the protomagnetar model (e.g., [6]) and the collapsar sce- 
nario [7] [8]. In the latter, the CCSN fails and a black hole 
(BH) with an accretion disk forms or a weak explosion 
occurs leading to fallback and BH/disk formation. 

In this Letter, we address, for the first time in 3 + 1 
general relativity (GR), the formation of spinning BHs 
in failing CCSNe in the context of the collapsar sce- 
nario of LGRBs. Our full GR method allows us to self- 
consistently follow core collapse, bounce, postbounce evo- 
lution, protoneutron star (PNS) collapse, BH formation, 
and the subsequent early hyperaccretion phase. For the 
first time, we extract the gravitational wave (GW) signa- 
ture of a failing CCSN that evolves into a collapsar and 
track the properties of the nascent BH with the dynam- 
ical horizon formalism [9]. 

Previous work on BH formation in CCSN/LGRB pro- 
genitors was limited to spherical symmetry [T0llT2] and, 
due to gauge choices, simulations could not be continued 
beyond BH formation. Multi-D studies either consid- 
ered isolated NS collapse (e.g., [13]) or BH formation in 
very massive polytropes [14] [15] . Recently, Sekiguchi & 
Shibata [16] carried out the first axisymmetric (2D) GR 
simulation that continued beyond BH formation in a hot 
polytrope, but did not extract the GW signal. 

Method. We employ the Zelmani 3 + 1 GR core col- 



lapse simulation package [17] which is based on the Cac- 
tus framework and the Carpet adaptive mesh refinement 
(AMR) driver [18], and uses the open-source Einstein- 
Toolkit for GR curvature (via [19]) and hydrodynam- 
ics evolution (via an updated variant of [13 ). We ex- 
tract GWs directly from the spacetime fields using the 
fully gauge- invariant Cauchy- Characteristic Extraction 
method of [I7j [20]. The simulations are performed in 
an octant of the Cartesian 3D cube with periodic bound- 
aries on two of the inner faces of the octant and reflec- 
tion symmetry about the equatorial plane. This limits 
3D structure to even £ and m that are multiples of 4. We 
use 11 levels of AMR, adding levels during collapse and 
postbounce evolution when needed. In our baseline reso- 
lution (BR), the finest resolution is ~370 m and ~92 m at 
bounce and BH formation, respectively. We also perform 
calculations with 20% higher/lower (HR/LR) resolutions 
and check stability and consistency by monitoring the 
ADM constraints. They show 2 nd -order convergence up 
to bounce and l st -order afterwards. After BH formation, 
convergence is reduced near the singularity, but the sim- 
ulations remain consistent and stable. ADM mass and 
angular momentum are conserved to < 3% in BR runs. 
All runs are carried out past BH formation, but only LR 
runs are continued to tens of ms after BH formation. 

We employ a hybrid polytropic-r-law equation of state 
(EOS; e.g., [21]). It smoothly matches a polytrope de- 
scribed by Ti « 4/3 at subnuclear densities with one de- 
scribed by T 2 > Ti at supernuclear densities, allowing to 
capture the stiffening of the nuclear EOS. A T-law com- 
ponent (described by T t h) accounts for thermal pressure 
contributions due to shock heating. We set Ti = 1.31 in 
the collapse phase and choose a rather soft supernuclear 
EOS by setting r 2 = 2.4. This results in a maximum 
non-spinning PNS gravitational mass of ~ 1.7 Mq, which 
provides for rapid BH formation, but is below the empir- 
ical NS mass limit [22]. We choose T th =4/3 for the 
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FIG. 1: Progenitor model u75. Left ordinate: Radial den- 
sity distribution. Right ordinate: Angular velocity (red dash- 
dotted curve) and specific angular momentum (blue curve) 
as a function of cylindrical radius w as given by our rotation 
law, reproducing features seen in the rotating progenitors of 
[4]. Vertical lines mark the enclosed mass. 
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postshock flow whose effective T is reduced by the disso- 
ciation of Fe-group nuclei. Neutrino heating (unlikely to 
be dynamically relevant in this scenario) is neglected, but 
we account for postbounce neutrino cooling of the outer 
PNS and the postshock region via the cooling function 
given in [23 . 

Initial Conditions. We use the 75-M , 10 _4 -solar 
metallicity model u75 of [24 whose compact core fa- 
vors early BH formation [12] . u75 could be a viable 
GRB progenitor if mass transfer to a binary compan- 
ion removed its H/He envelopes. We map u75's inner 
~5700km (enclosed mass ~4.5M ) onto our 3D grid 
and impose constant rotation on cylindrical shells with 
radius w via a rotation law motivated by the GRB pro- 
genitors of [4]: The inner iron core is in near uniform 
rotation and ft drops oc w~ 2 further out. Close to the 
edge of the iron core, Vt drops by a factor of order unity, 
then continues to decrease oc r -( % with < £ < 2, 
leading to a radial increase in the specific angular mo- 
mentum j, endowing mantle material with sufficient spin 
to form a disk at small radii. The functional form is 

n(w) = (i - \(w))n (i + (w/A) 2 )- 1 + \(w)^q(i + 

(zut/A) 2 )- 1 ^ + (max(0,tu - m t )/A))S. Here, \{w) = 
(1 + tanh((tu - w t )/5w))/2. We set A = 1000km, 
wt = 1950 km, £ = 1/3, and 8w = 100 km. Qq is the 
central angular velocity that we vary from to 2rads _1 . 
Fig. [I] depicts u75's density profile along with Q{w) and 
for the = 2rads _1 case. Model names, param- 
eters, and key results are given in Table [I] 

Dynamics. — The homologous collapse of the inner core 
to nuclear densities proceeds as in the standard CCSN 
case. For the initial inner core rotation rates considered 
here, centrifugal effects are negligible in the prebounce 
phase and all models reach core bounce after ~ 114 ms of 
collapse. A hydrodynamic bounce shock is launched, but, 
due to neutrino cooling and the low V in the postshock 



FIG. 2: Top: Maximum density y9 ma x and central ADM lapse 
function a m in as a function of postbounce time in all models. 
After horizon formation, the region interior to it is excluded 
from min/max finding. Bottom: BH mass and dimensionless 
spin a* as a function of postbounce time. All models follow 
the same accretion history once a BH forms and settles down. 
The data shown in this figure are from the LR runs since these 
were carried out longest after BH formation. 



TABLE I: Model summary. Qo is the initial central angu- 
lar velocity. £bh is the time after bounce to BH formation. 
M max is the mass of the PNS at that time, a* and a* are the 
dimensionless BH spin shortly after BH formation and when 
the simulation is stopped, respectively. Eqw is the emitted 
GW energy and f c is the characteristic GW frequency [25] in 
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region, quickly (within milliseconds) succumbs to the ram 
pressure of the outer core, which is accreting at a rate of 
initially tens of M s _1 . The shock stalls at only <50km 
and gradually retracts in all models. In the top panel of 
Fig. |2j we plot the maximum rest mass density p max (t) 
that rapidly increases as accreted material settles onto 
the outer PNS core. The slope of p max is steepest in the 
nonrotating model whose PNS becomes unstable earliest. 
In rotating models, centrifugal effects lead to an oblate 
and less compact PNS that contracts more slowly and is 
stable to larger mass (cf. [12] [26] and Table 0. The time 
to BH formation and the maximum PNS mass increase 
roughly with Qq. 

Once dynamical PNS collapse sets in, an apparent 
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FIG. 3: Snapshots of the meridional density distribution with 
superposed velocity vectors in model u75rotl taken at various 
times. The top left panel (note its special spatial range) shows 
a snapshot from 10 ms after bounce. The top right and bot- 
tom left panels show the point of PNS instability and the time 
at which the AH first appears, respectively. The bottom right 
panel, generated with a separate color range, shows the hy- 
peraccreting BH at ~ 15 ms after its formation. All colormaps 
have density isocontours superposed at densities (from outer 
to inner) of p = (0.1, 0.25, 0.5, 0.75, 1.0, 2.5, 5.0) x 10 10 gem -3 . 



horizon (AH) appears within —lms and quickly engulfs 
the entire PNS. With the PNS and pressure support re- 
moved, postshock material and the shock itself immedi- 
ately subside into the nascent BH. The bottom panel of 
Fig. [2] shows the evolution of BH mass and dimensionless 
spin a* in all models. The former jumps up as the AH 
swallows the PNS and postshock region, then increases 
at the rate of accretion set by progenitor structure and 
is largely unaffected by rotation at early times. The di- 
mensionless spin reaches a local maximum when the BH 
has swallowed the PNS core, then rapidly decreases as 
surrounding lower- j material plunges into the BH. This 
is a consequence of the drop of j at a mass coordinate 
close to the initial BH mass (cf. Fig.|TJ. Table summa- 
rizes for all models the values of a* at its peak and at the 
time we stop the LR run. 

In Fig.[3j we plot colormaps of the density in the merid- 
ional plane of the spinning model u75rotl taken at var- 
ious postbounce times. The rotational flattening of the 
PNS is significant and so is the centrifugal double-lobed 
structure of the post-BH-formation hyperaccretion flow. 
The latter is unshocked and far sub-Keplerian with in- 
flow speeds of up to 0.5c near the horizon. The flow will 
be shocked again only when material with sufficiently 
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FIG. 4: Top: GW signals /i+ >e emitted by the rotating models 
as seen by an equatorial observer and rescaled by observer 
distance D. Bottom: Spectrogram of the GW signal emitted 
by the most rapidly spinning model u75rot2. 



high specific angular momentum to be partly or fully cen- 
trifugally supported reaches small radii (cf. [H]). Based 
on progenitor structure, our choice of rotation law, and 
the assumption of near free fall, we estimate that this 
will occur after —1.4 s, —2.4 s, —3.9 s in model u75rot2, 
u75rotl.5, u75rotl, respectively. At these times, the 
BHs, in the same order, will have a mass (a*) of — 8M 
(0.75), ~14M© (0.73), and -23 M (0.62). 

GW Signature. — The top panel of Fig. [4] depicts the 
GW signals emitted by our rotating models. Due to the 
assumed octant symmetry, GW emission occurs in the 
I = 2,m = mode. The nonrotating model leads to 
a very weak GW signal and is excluded. At bounce, a 
strong burst of GWs is emitted with the typical signal 
morphology of rotating core collapse (e.g., [27 ) and the 
peak amplitude is roughly proportional to model spin. 
Once the bounce burst has ebbed, the signal is domi- 
nated by emission from turbulence behind the shock. It 
is driven first by the negative entropy gradient left by the 
stalling shock and then by neutrino cooling, whose effect 
may be overestimated by our simple treatment. Interest- 
ingly, the signal strength increases with spin. This is not 
expected in a rapidly spinning ordinary 2D CCSN, since 
a positive j gradient in the extended postshock region 
stabilizes convection. In our models, the postshock re- 
gion is considerably smaller and shrinks with postbounce 
time. The driving entropy gradients are steeper and the 
change of j in the postshock region is smaller. Also, in 
contrast to 2D, our 3D models allow high- mode nonax- 
isymmetric circulation. We surmise that the combination 
of these features with increasing spin (feeding greater cir- 
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culation) results in a stronger GW signal. 

The intermittent period of turbulent, low- amplitude 
GW emission ends when PNS collapse sets in, leading 
to a second pronounced spike in the waveform, marking 
BH formation. The collapse signal evolves into the ring- 
down emission of the nascent BH that rapidly assumes 
Kerr shape. The GW emission ceases soon after and the 
un-shocked axisymmetric accretion flow does not excite, 
at appreciable amplitude, BH quasi-normal modes that 
could emit GWs. The strength of the BH formation sig- 
nal scales with and its gLE^gw / df peaks at ~3.9 kHz, 
-3.4 kHz, -2.9 kHz, in u75rotl, u75rotl.5, and u75rot2, 
respectively. The lower panel of Fig. [4] shows the spectro- 
gram of the GW signal in model u75rot2. There is a clear 
trend towards higher frequencies during the postbounce 
pre-BH phase, but BH formation itself, while peaking 
in the kHz range, leads to significant emission also at 
lower frequencies, which is favorable for detection by ad- 
vanced laser-interferometer GW observatories (aLIGOs). 
In Table [T| we provide quantitative results on the GW 
emission in our model set. For an event at lOkpc, we 
estimate optimal single-detector aLIGO signal-to-noise 
ratios (see [HI [25]) of -36 (u75rotl), -68 (u75rotl.5), 
and —94 (u75rot2), and —6 for the nonrotating model 
u75rot0. Note that real GW burst searches will not re- 
cover all available signal power. 

Discussion. — We have performed self-consistent 3 + 1 
GR simulations of stellar collapse in the context of the 
collapsar scenario for LGRBs. Albeit approximate in 
many aspects, our models elucidate characteristic qual- 
itative features in the dynamics and GW signature of 
these events. The rotating-collapse-bounce-PNS-phase- 
BH-formation-hyperaccretion sequence and its GW sig- 
nature are robust aspects of the early collapsar evolution. 
More realistic physics will undoubtedly affect quantita- 
tive results, but the overall qualitative picture is unlikely 
to change. The characteristic GW signature seen in our 
models will enable aLIGO to distinguish between a suc- 
cessful and failed CCSN purely on the basis of observed 
GWs, provided the event is sufficiently nearby. 

A more realistic, stiffer EOS will increase the delay be- 
tween bounce and BH formation and will lead to higher- 
amplitude, lower-frequency GWs. An improved neutrino 
treatment may reduce the vigor of turbulence in the PNS 
phase and decrease the amplitudes of the associated GW 
signal. Symmetry- free 3D evolution could reveal nonax- 
isymmetric dynamics that may lead to an enhanced GW 
signal [27]. Only the inclusion of MHD may lead to a 
large qualitative change by potentially leading to a strong 
explosion, leaving behind a magnetar [6j [28j [29]. This 
study is a first pioneering step and much work lies ahead 
before a clear and quantitative picture of the CCSN- 
LGRB connection can be drawn. 
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